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1 Introduction 



Reaction-diffusion processes are the subject of much research |jj || [25] [|1| 
H [23], a reaction-diffusion process occurs as reactants in a solution dif- 
fuse in the liquid and react amongst themselves. A common approach to 
reaction-diffusion processes is to consider the density fields of the different 
reactants participating in the reactions. This approach stands in contrast 
to the more naive approach of tracking the locations of the different reac- 
tants, or computing the wave functions of the different reactants. Whatever 
approach is taken the interest in a reaction-diffusion system is usually in its 
spatio-temporal evolution. The density field approach is especially adept 
for this purpose, since the actual location of specific reactants is, usually of 
no interest. In the density fields approach the spatio-temporal evolution is 
modeled through partial differential equations (PDE's). 

Another approach to reaction-diffusion processes that we have suggested 
is the microscopic approach. In this approach we consider the number of 
reactants at discrete lattice points, where the lattice models space. The 
main difference from the density field approach is that rather than using 
continuous densities in a continuous space as in the density field approach, 
we use discrete densities in discrete space. 

The microscopic simulation approach is closer to the real simulated sys- 
tem when there are only trace densities of the different reactants. This is 
because it is in this situation that the discrete nature of the reactants comes 
into play. Consequently the PDE approach describes the system with less 
accuracy than when there are many reactants. 

Reaction-Diffusion processes are not restricted to describing chemical 
systems. Indeed reaction-diffusion processes have even been used extensively 
in population biology |2l}| . We have also used a reaction-diffusion model in 
a marketing context. We have seen that discretization was crucial in the 
behavior of the modeled market. Thus showing that microscopic simulation 
could be of use to researchers who need to model real-life systems. 



2 The Density Field Approach 

2.1 Analytical approach 



In this section we shall describe the density field approach [21]. As men- 
tioned in the introduction, The density field approach considers the evolu- 
tion of the the density fields of the reactants participating in the reaction- 
diffusion system. Say that reactants numbered 1 to n are participating in 
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the reactions (possibly as reactants or as products). Then we could label 
the density fields by pi(x,t). The density is defined as: 

f ^.s 1 r N(i,A,x,t) 

Pi(x,t) = — hm — — rr (1) 

n v(A)->o V(A) 

A is a box, N(l, A,x,t) is the number of molecules of type i that are 
in the box A located about x. V{A) is the volume of box A. And no is a 
constant that serves the same purpose as Avogadro's constant. It should 
be noted that, if a smooth density function is wanted, the limit should not 
be taken to zero literally, but rather should be taken down to a scale much 
larger than one that shows the discretization of the reactants, and much 
smaller than the scale of the macroscopic spatial patterns. 

Let us assume I possible reactions, where reaction i is of the form: 

S h . + ... + Sj . — ► S k , - + ... + S k .. (2) 

JX,l Jmj,i ^l,z "''l*,? V / 

Where S r denotes the r th reactant. As an example of such a reaction let us 
look at: 

S 1 + S 1 +S 2 ^S 1 + S 1 + S 1 . (3) 

this reaction is of the form (||) . An interpretation of this reaction is that 
two reactants of species 1 can cause a reactant of species 2 to turn into a 
reactant of species 1. Let us consider a system which has two chemicals, 
which we shall denote, as usual, by S± and S 2 - These chemicals can diffuse 
with diffusion coefficients of D\ and D 2 respectively. These chemicals can 
also react according to the reaction scheme (||). 

The time evolution of the fields, pi(x,t) and p 2 (x,t), is given by: 

^ = Dl V 2 P i + k ■ p\ P2 (4) 

?§ = D 2 V 2 p 2 -k-pjp 2 (5) 

Equation |] has two terms on the right hand side (RHS), let us turn our 
attention first to the second term. The term contains the expression p\p 2 . 
This is proportional to the probability that two reactants of species 1 and 
one reactant of species 2 meet in a small region in space (the volume of 
that region is a given). The coefficient k is the probability that, once the 
reactants met in the small region in space, they will react with one another. 
The Diffusion term is the familiar term, which originates from the "random 
walk" of the reactants. 
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In general, when we have I reactions, each of the form (^), 



+ ••• + Sj, 



Sk lti + ■•• + Sk, 



(6) 



If we denote by N Pir the number of reactants of species r created or 
annihilated by reaction p then we have the following rate equations: 



2.2 Simulation by finite difference 

The former section introduced the formalism of the density field approach 
which is essentially analytic, but there is no data structure on a computer 
that can hold an arbitrary continuous field. So space is discretized in the 
computer simulation. The other problem which arises is the need to integrate 
the differential equations over time. This is done again by discretization but 
the solution now is to discretize time. The most naive way to integrate using 
the differential equation is by Euler integration. This method's main draw- 
back is the computation time that it requires to get accurate results. But 
fundamentally it is no different than other more sophisticated methods such 
as the runga-cutta method. We shall outline this finite difference approach 
using Euler integration below. 

In the finite difference approach using Euler integration one replaces the 
fields p r {x, t) with x £ R d and t E R by p*(x, t) where x £ uj d t G uj (to being 
the natural numbers). Let us suppose that the following equations hold for 
the fields p r : 



In order to make this transition from p to p*, space is conceptually 
divided to a discrete d-dimensional lattice of spacing Ax and time is divided 
to a discrete 1-dimensional lattice (actually a series) of spacing At. Now the 
aim is to make the following equality be a good approximation: 



The way to make this approximation good is to let Ax and At be small 
and to let p* follow the dynamics: 







p*((n 1 ,n 2 ,...,n d ),m) « p r (Ax • (nx,n 2 , ...,n d ) ,m ■ At) 



0) 



p* r (x,t+l)-p* r (x,t)=At-f(pt,...,p* n ) + -^- 
(p* r (x + (1, 0, ...0), t) + p* r (x + (-1,0, ...0), t) + 
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p* r (x + (0,1,. ..0),t)+p* r (x + (0,-1,.. .0),t)+.. . 
p* r {x + (0, 0, ...1), t) + p* r {x + (0, 0, ... - 1), t) - 2dp* r {x, t)) (10) 



The first term on the right hand side of the equation is the first order 
approximation of the difference between p r {x,t) and p r (x,t + At), after a 
time interval of At has passed assuming the dynamics (^) . The second term 
accounts for diffusion and includes the discretization of the V 2 operator. 
This term includes positive and negative terms. The positive terms are 
contributions to the density at site x from densities at neighboring sites. 
Neighboring sites are those sites which have the same coordinates as x but 
for a single coordinate which must be only one lattice point away. This 
contribution is due the fact that diffusion causes chemicals to move from one 
location to another in space. The negative term accounts for the chemicals 
leaving site x. 

Now, if we replace / with the terms from the density field approach to 
reaction-diffusion, we get: 

p*(x,t + 1) - p*(x,t) = 

l s=m p 

At^VV 1] p* s Jx,t) 

p=l s=l 

+^^-{ P * r {£ + (!> °» -o)»*) + + (-!> °» -o)»<) + 

p* r (x+(Q,l,...0),t)+p* r (x + (0,-l,...0),t)+... 

P ;(x + (o,o,...i),t)+ P ;(x + (o,o,...-i),t)-2d P ;(x,t)) (ii) 



3 Microscopic Simulation of Reaction-Processes 
3.1 Fundamentals of the Approach 

In the previous section we have seen the prevailing approach for dealing 
with reaction-diffusion processes. Another approach that can be used is the 
microscopic simulation approach. This approach takes discretization one 
step further in the sense that the fields are discretized themselves, but takes 
a welcomed step backwards in the sense that time is not discretized. This 
approach is useful because it is closer to reality. Chemicals are discrete 
entities (at least in the classical approach which is a good approximation in 
solutions). 

Again we have a field p**(x,t) where x G uj d t G lo, but this time 
p**(x,t) G lo. We have said that time is not considered discrete in the 
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microscopic simulation model, but nevertheless we have t £ u. This is not 
a contradiction it is simply an expression of the fact that reactions occur 
at discrete time points. Let us denote simulation discrete time with t* and 
real time with t. Say the reactions occur, in the real system, at times t n . 
Then when the simulation is at time t* = n it is supposed to approximate 
the real system at time t n . Simulating the real system between the times t n 
is of no use since nothing happens, except for diffusion which we should also 
treat as a reaction. But the problem is that diffusion, in the real system, is 
not a process that occurs at time points, rather it is a continuous process. 
On the other hand we have modeled space by discrete sites. Diffusion is 
modeled by chemicals hopping from one site to another. This process is 
discrete since chemicals hop at discrete time points. So amongst the times 
t n there are times at which the reaction which takes place is diffusion, that 
is to say hopping of chemicals to neighboring sites. We should stress that 
the time interval between t n and t n+ \ is not a constant. So the real time is 
not approximated by t* ■ At for some At. 

The time interval between t n and i n +i is large when the time interval 
between successive reactions, in the real system, is large. Roughly speaking 
this happens when there are not many reactants, or many inert reactants. 
This is also when the microscopic simulation is at its best (in terms of the 
simulation's speed), since the simulation's single step covers a lot of time. 
We shall give a quantative estimate for this time interval later. 

Let us now turn to the relation of p** to the real system that it is 
supposed to approximate. We assume space of dimensions d. Let us denote 
by N r (x, I, t) the number of reactants, in the real system, of species r, located 
in a box of length / around x at time t. The approximation relation is given 
below: 

p**(x,t*) ~ N r (Ax-x,Ax,tt*) ^p* r (x,t t *) ■ Ax d n (12) 

Let us imagine a grid of spacing Ax dissecting the real system, so that 
space is divided into little boxes. This division is not physical, but mental. 
Now each such box is simulated as a site on the simulation lattice. The 
number of reactants at each lattice point should approximate the number of 
reactants in a the little boxes imagined in the real system. This is the nature 
of the first approximate equality in equation (]i~2l ) . The second approximate 



equality in equation (12) is due to the approximation of the finite difference 



approach to the real system. 
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3.2 The Monte-Carlo Method 



Now we shall introduce the dynamics of p**. Since the microscopic simula- 
tion is a Monte-Carlo simulation, p**'s dynamics follow the following rules: 

1. Choose p**(x,t* + 1) from a probability space, f2o(i*). 

2. Calculate the new probability space Qo(t* + 1). 

f2o(i) associates a probability for every possible p**(x, t+1), but actually 
we are intent on performing one reaction at a time. So ilo(t) will give non- 
zero probabilities only for those p**(x,t + 1) that differ from p**(x,t) by a 
single reaction (or diffusive hopping). So we can look at ilo(t) as associating 
a probability for every possible reaction. So let us construct a probability 
space fi(9fc, t) which associates a probability for each possible reaction. 

Let us speak of a reaction, -ft. This reaction can potentially take place 
anywhere in real space. In particular the reaction can fall within any one 
of the little boxes that we have discussed in the previous sub-section. We 
shall denote reaction 3ft taking place at a little box corresponding to site x 
by Dftj. Our task now is to find out, given some initial conditions, what is 
the probability that the reaction that will take place next is Dftj. 

Let us expand a bit on the stochastic process that the real system under- 
goes. The stochastic process is comprised of events (reaction and diffusion) 
occurring stochastically at discrete time points. The events we are consider- 
ing are the real reactions and the movement of reactants from one little box 
to an adjacent one. Let dt be a small time interval, then for 9ft^ there is a 
chance P($lg,t*)dt that this reaction will occur in the time interval dt. The 
probability density, P($tg,t*), is called the reaction rate for reaction 3ft^, 
and it is exactly what we used in order to formulate the PDE for the real 
system, as we shall soon see. First let us notice that this probability density 
has reciprocal time as units, which is consistent with the name 'rate'. 

The connection of P($tg,t*)dt to the PDE's will be useful in calculating 
PQR-S, t*). So let us explore this connection. Reaction Oft is of the form (1), 
that is to say: 

+ ••• + Sjmx,x *■ ^ fc i,» + ••• + ^n s ,»- ( 13 ) 

A crucial assumption for the validity of the PDE's is that there is some 
time scale, dt, during which P($tg, t*) does not change much and still for 
the same time scale, dt, many reactions occur. Under these assumptions it 
can be shown that the number of reaction and diffusion events that occur 
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at the time interval dt in the box around x is, simply, P(^,t*)dt. Let us 
look at species r ( recall the notation used in equation (|^) ). Now letQ 

Nn,T = #{n\n < n^,k n ^ = r} - #{n\n < m^j n ^ = r}. (14) 

This NsR r gives the number of reactants of species r created in the reac- 
tion 3? ( a negative number indicates that the species is annihilated in the 
reaction). So the number of reactants of species r formed, at the box around 
x, in the time interval dt is: 

P(^,t*)-dt-N^ r (15) 

The number of reactants of species r at the box which corresponds to 
site x at time t n is approximated by p**(x, n). So the rate at which p**(x, n) 
changes due to reaction 3?^ is given by: 

=P(dt s ,t*)-Nn, r . (16) 

Where the index 3J^ denotes that the reference is to the rate of change 
due to reaction 31g alone. We have already seen the rate of change in the 
finite difference case (equation (|TT1)). There we expressed as a sum of 
terms, each expressing a rate due to different reactions. Another term was 
due to diffusion. Realizing that the terms appearing in the finite difference 
case express the same thing as the rate expressed at ( |i~6| ) modulo the approx- 



imation relation (12) we can find an expression for P(3?£, t), this is given 
by: 



P(fe t f) = Ax d n ■ % J] (Ca^J 

This assumes that the reaction in question is a real reaction, that is 
not diffusion. For the case of diffusion of species I we have the following 
equation: 

P(ft s ,t*) = 2d-D r j¥- S (18) 

This equation is derived considering the last term in equation (|IT|), 
— 2d fV At p* (x,t), which is due to reactants hopping from site x to neigh- 
boring sites. And then considering that P(dtg,t*), in equation (|l8|), is the 
probability density for hopping to neighboring sites. After using the approx- 



imation relation (12) we get (FT 



1 #S denotes the number of elements in S. 
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Now that we have an expression for the rates P(3ig, t*) of the different 
reactions, we still face the task of finding the probability associated with each 
possible reaction (including diffusion) by f2. Again let us point out that the 
probability that f2 associates with reaction at time t* is the probability 
that this reaction will come next in the sequence of reactions. Now between 
the times tf* and tt*+i nothing happens in the reaction chamber, apart from 
reactants moving inside the little boxes we have imagined. In this time 
interval the reactants don't cross the boundaries of the boxes. If the rate of 
is twice that of ift^ then we should expect that reaction has twice 
the chance to be the next reaction that occurs than 0?^2- Let us expand 
a bit on the nature of the assumption we made in the previous statement. 
We assume that the probability distribution associated with Q is dependent 
only on the situation of the current configuration of the reaction chamber 
and is independent of the time that passed since the last reaction-diffusion 
event. So it is actually this assumption (that the stochastic process has no 
memory) that allows us to compute the probability of different reaction and 
diffusion events associated by Q as a function of the rates. The preceding 
argument leads to the following relation, given two reaction-diffusion events 

i and 2 : 

m^t*) p(sfe )2 ,t*) 1 j 

Taking into account the normalization of probability distribution to one, 
we now can compute the probability associated by ft to a reaction-diffusion 
event $tg,t*. The solution is: 

m^n = T m P % n t * y (20) 

4 Anderson Localization in a React ion- Diffusion 
System 

4.1 The Reaction-Diffusion System 

Anderson localization is a phenomenon associated with electron-transport 
behavior in disordered materials. The disorder in the material induces a 
phase transition of the electron eigen-functions from an unlocalized state in 
which Ohm's law is valid into a localized state in which the material behaves 
as an insulator. The wave functions of an electron under the influence of a 
periodic potential is periodic. A question arises, in the context of disordered 
materials, of how this periodic wave function is influenced by disturbances 
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to the periodicity of the potential. One might think that the eigenfunction 
of a slightly perturbed potential (that is to say perturbed from an origi- 
nally periodic potential), would be eigen-fucntions slightly perturbed from 
a periodic function. This intuition proves misleading in the metal-insulator 
transition case. It is seen that some of the eigen-functions exhibit a marked 
departure from periodic functions, even for small disturbance of the poten- 
tial. These functions are seen to be localized. Meaning that there are small 
"islands" in which the wave function has high modulus and there are large 
spaces between these islands where the wave-function has low modulus. In- 
deed the wave functions that depart from periodicity have the approximate 
form : 

e~ « (21) 

£ is the localization length of the eigenfunction. This approximation is 
good for the tails of the eigenfunction. 

This Anderson localization effect is also observed in reaction diffusion 
system. Let us take the example given by Shnerb and Nelson. They pre- 
sented a reaction-diffusion system described the schematics: 

A + F — > A + A + F 

A + A — >0 (22) 

A is the only reactant to undergo diffusion. 
Where signifies that the reaction has no products (alternatively it can be 
interpreted that the products of the reaction are inert). These reactions to- 
gether with diffusion can be seen as the schematics of the population biology 
of bacteria. The first reaction accounts for the reproduction of the bacteria. 
The rate at which the bacteria reproduce is controlled also by the concen- 
tration of F. F can represent, for example, the intensity of light that falls on 
the bacteria pQ ]. But can also represent any factor that controls the rate 
of bacteria reproduction, with the provision that this factor is not variable 
in time, and in particular non-exhaustible. The second reaction accounts 
for the dying of bacteria due to overcrowding. The reaction dynamics has 
the species A facilitating the production of more A's. This property A is 
called autocatalysis. Autocatalysis is an important concept in pattern for- 
mation f23|l , the origins of life || [18] Q [10] Q and economy |22| (where 



autoctalysis is the underlying concept of multiplicative dynamics) . 

To see the similarity of this problem to an eigenfunction problem in 
quantum mechanics, such as Anderson localization, let us write down the 
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PDE's which are associated with the reaction dynamics above: 

D A V 2 PA + PA- PF- Pa (23) 



dt 

which can be reformulated into 
dpA 



nt (D A V 2 + PF - Pa )(pa) (24) 
Now if we drop the last term, we have: 

^ = (D A V 2 + p F ){p A ) (25) 
which is quite analogous to the time dependent Schrddinger equation: 

{ W -<5^ + <7(*.«»M (26) 

Where we have U playing the part of pp, ^ pl a yi n g the part of D A and 
ip playing the part of p A . We have only to remember that we have dropped 
the quadratic term and that there is an additional coefficient, -i, in the 
Schrddinger equation 0. Consequently, the reaction-diffusion PDE, without 
the quadratic term, can be seen as a Schrddinger equation with imaginary 
time. 

Notice that we did not write an equation for since pf does not 
undergo any dynamics and therefore does not change in time. The phe- 
nomena of Anderson localization determines that if pp is not constant in 
space, but rather is stochastic, then we shall have these localized states that 
we have described above, both for the quantum mechanics case and for the 
reaction-diffusion case. In the quantum mechanics case, a constant or peri- 
odic potential, U(x,t), entails a periodic field, if). In the reaction-diffusion 
case, a constant potential, pp, entails a constant field, p A . 

A constant, stochastic potential naturally arises when speaking of disor- 



dered materials [19]. But in reaction-diffusion systems, such a potential is 
more of a stretch. We have given the example of light intensity as a stochas- 
tic constant potential. But in many more cases the reaction rate associated 
with the reproduction of bacteria, or some chemical, is dynamic. This leads 
us to deal with a dynamic potential, U(x,t). Nelson and Shnerb have al- 
ready discussed a time-dependent potential of the form U(x — vt) (actually 
they have looked at the case where the media is moving, which is equivalent 
to a moving potential in the opposite direction). 

2 i here is 



11 



4.2 Results of microscopic simulations for Anderson local- 
ization 

We simulated (using microscopic simulation) the reaction dynamics associ- 
ated with Anderson localization . The reaction schematics (together with 
the reaction rates) is given belowQ: 

1. A + C^ A + A + C(15) 

2. A + A^0(3) 

3. A -> 0(60) 

The size of the reaction chamber is 300x300 lattice points. The average 
number of molecules of type C per site is 5. They are dispersed in the 
beginning of the simulation randomly. That is to say that the positions 
of the 5x300x300 molecules of type C are chosen at random. The spatial 
distribution that resulted can be seen in figure ([!]). The initial distribution 
of A is similarly dispersed but with and average of 2 molecules of type A 
per site. Only A molecules diffuse. The simulation's diffusion rate is 3. 
After awhile a steady state is reached. Off-course C is not dynamic, and it 
remains as in figure (|l]). On the other hand A is dynamic and the steady 
state is depicted in figure ((D) 

4.3 Anderson Localization - in the search of dynamic clus- 
tering 



Shnerb and Nelson [2C] have explored the consequences of changing the po- 
tential U (x) by a moving potential U(x — vt), what would happen if instead 
of having the potential drift at a constant speed the potential would undergo 
diffusion itself? This would correspond to replacing ( p3| ) by the two coupled 
equations: 

d PA _ r, . w2 2 



A4V PA + PA - PF - Pa 

D F V 2 PF (27) 



dt 

dp F 



Of, 

We have found no treatment in the literature for this kind of moving 
potential for a good reason: if pp is treated as non-negative continuous 



The numbers in parentheses at the right of the reactions denote their corresponding 
rates 
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Figure 1: Snapshot of molecules of type C. 
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Figure 2: Snapshot of molecules of type A. 
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variable as is usually done in the context of differential equations and since 
the only dynamics imposed on pp by (|27|) is diffusion, pp converges to a 
steady spatially-uniform state pp = const, giving rise to the trivial solu- 
tions of a system with no potential fluctuations at all. But if the potential is 
treated as a discrete variable diffusion does not necessarily induce the uni- 
form distribution of the potential and localization effects still have a chance 
to prevail. This situation leads to clear differences between the different 
simulation approaches we have described, in further sections we will show 
other examples in which discretization leads to different results in different 
simulation methods. 

Extending the analogy to population biology we could look at A as repre- 
senting a population of parasites dependent on a host species F for a-sexual 
reproduction. The population represented by F is not affected by the para- 
sites and performs random diffusion in space. The phenomena we are most 
interested in is dynamic clustering or grazing. Dynamic clustering would 
under our analogy represent parasite herds moving in space due to changes 
in the spatial distribution of the species they need in order to reproduce. 

Translating our situation into a reaction-diffusion system will give the 
same results as (p^) but in this case both A and F undergo diffusion. We 
have added to (|22|) the reaction: 

A — ► (28) 

representing the dying of A not due to overcrowding. 

Since ( [22] ) describes no creation or elimination of F the total number of 
Fs in all the lattice sites is constant throughout the simulation and < F > 
- the average concentration of F is constant as well. Given a constant value 
for Da the results of simulations of (|29|) are controlled by < F > and Dp 
We have simulated (using microscopic simulation) the following reaction- 
diffusion system with different values of < F > and Dp and keeping Da = 6: 

B 

A + F -> A + A + F(13) 
A + A -> 0(0.01) 

A -» 0(8) 

(29) 

The size of the simulation is 128x128. The system was seeded with 3 
reactants of species F per cell and 1 reactant of species A per cell (the 
reactants were placed at random locations). 
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Figure 3: Snapshot of A's concentration, simulation parameters are Da = 
6,D F = 4,< F >=0.15) . 




Figure 4: Snapshot of the same simulation as in (|3|) at a later time, the 
clusters have moved. The result of this simulation falls into the (3 category 
in our notation. 
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Figure 5: A's fill up the simulation space, this corresponds to a 7 situation 
in our notation. 

We have divided the results of the simulation into 3 categories: 

• a- The simulation results in A filling the whole simulation space. 

• (3- The simulation results in A forming dynamic clusters. 

• 7-The simulation results in the total extinction of A from the simula- 
tion space. 

In fig(^) we can see an abrupt phase-transition between a states to 7 
states for high Dp values around < F >= 0.61. We shall give a theoreti- 
cal explanation to this phase change. For high Dp values F reactants take 
shorter times to pass between different parts of the simulation space. There- 
fore we can assume that A reactants are affected by all F reactants in the 
simulation space and the mean-field approximation: 

F=<F> (30) 

is valid. 

4 Numbers in parentheses at the right of reactions are the corresponding reaction-rates. 
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Figure 6: Results of simulations of the system(^) for different values of 
< F > (x-axis) and Dp (y-axis) . Notice that dynamic clustering occurs also 
in the realistic range of: 0.5 Da < Dp < 1.5-Da- Circles denote simula- 
tion resulting in a situations, Asterisks denote simulations resulting in (3 
situations and squares denote simulations resulting in 7 situations. 



18 



302 




238 



0.50 



0.55 



0.60 



0.65 



0.70 



Figure 7: Results of simulations of the system(p9|) for Dp = 300 (y-axis) 
as a function of < F >(x-axis), circles denote simulations ending in an a 
situation and squares denote simulations ending in a 7 situation. The change 
in behavior is abrupt and occurs in the 0.55-0.61 region. Simulations in that 
parameter region resulted in (5 situations. 
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Assuming (30) we are now looking for solutions to: 



^ = D A V 2 p A + k lPA • < F > -k 2 p\ - k 3PA (31) 
Keeping in mind that A can take only integral values and dropping the 



diffusion term from (31) we find that: 

&2 + h 



A = 1 (32) 



< F >-- 



is an unstable node fixed point of (|3l|). This leads to predict the death 
of As for < F > smaller than k2 ^ 3 and that As fill up the space for < F > 
larger than fc2 + fc3 . 

In the case of the simulated set of reactions (^) we have: 

« 0.61 (33) 

h 

The preceding argument explains the sharp transition between a and 7 
states for the high Dp = 300 value as can be seen in fig(0). 

One might be tempted to think that this mean-field argument will be 
sufficient in explaining the emergence of grazing at low Dp values: for low 
diffusion coefficients persistent spatial discrepancies in F's concentration 
prevail throughout the simulation space, some areas would be affected from 
a local F concentration larger than 0.61 and would sustain a population 
of As and some areas with a lower local F concentration would be empty 
of As. These spatial discrepancies change with time causing our clusters 
to move across the simulation space. If this was the sole mechanism lead- 
ing to the grazing behavior one would expect that the /3 simulation results 
would appear for low Dp values equally distributed on both sides of the 
< F >=0.61 asymptote. Clearly, fig(g) shows us that this is not the case, 
the (5 situations are concentrated around lower and lower < F > values as 
Dp decreases. Therefore we should search for another mechanism in order 
to explain the behavior seen in fig@. For lower Dp values, D A is not neg- 
ligible and clusters are not only supported by an influx of F reactants, but 
are also supported by their own ability to move and "find" areas of high F 
concentration in their surroundings. This mechanism of "searching" for F 
reactants should be part of an explanation for the behavior seen in fig(||) at 
low Dp values. 
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5 Microscopic simulation of the Gray-Scott model 



5.1 The Gray-Scott model - a pattern formation mechanism 

The Gray-Scott modelS was first designed as a model of glycolysis and as in 



the simplest form of Turing [23] pattern formation it involves two reactants 
one of them enhancing the auto-catalysis of the other, but the geometrical 
patterns resulting from this model are different from the ones observed in the 
Turing pattern case and unlike the Turing pattern case pattern formation 
occurs when diffusion coefficients are equal as well. The model is given by: 

^ = V 2 Pc -pcP 2 A + F(l-pc) , 

8pA ^ „2 2 



dt 



DaV 2 pa + PC Pa ~{F + k) PA (34) 



where pa (x, t) and pc (x, t) are the concentration fields for two chemical 
reactants A and C respectively. The terms F(l — pc) and — (F + k)pA in 
( |34|) describe the system as being in contact with an external reservoir kept 
at pc = 1 and pa = 0. Indeed (|34"|) accepts the solution: 



PC = 1 

PA = (35) 



Consider the equations that result from (34) by dropping the diffusion 
terms. The above mentioned stable fixed point exists throughout parameter 
space but for some range of F and k Pearson [15] has found another fixed 



point. Fixing k, increase or decrease of F causes the second steady state 
to be lost. The process of losing or gaining fixed-points as a function of 
the system's parameters is called bifurcation. Looking at the non-diffusive 
system's behavior under change in the external parameters can be very useful 
in order to understand the general behavior of the diffusive-system under 
microscopic simulation. 

Two types of bifurcations are usually distinguished - saddle-node bifur- 
cation and Hopf bifurcation. 

In a saddle-node bifurcation either a fixed point appears and splits into 
two fixed points or two fixed points become one and then disappear. The 
main feature of such bifurcations is the nature of the single fixed point at 
the bifurcation. The linearized system must have one zero eigenvalue and 
one non-zero eigenvalue at that point. In other words if 

^ = M/fi) (36) 
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is the linear expansion of the system around the fixed point, with M 
being the 2x2 coefficient matrix, then at the bifurcation point efei(M) = 
and tr(M) / 0. 

Hopf bifurcation occurs when an unstable (stable) focus goes through 
the creation of a limit cycle and becomes stable (unstable). At the bifurca- 
tion point both eigenvalues are purely imaginary and the Real part of the 
eigenvalues is positive (negative) before the bifurcation point and negative 
(positive) after the bifurcation point. 

In the case of our system, given k the second fixed point is lost through 
saddle-node bifurcation as F is increased and by Hopf bifurcation as F is 
decreased. 

Pearson [25] considers pc as the density of a liquid fuel and pa as a 



temperature field. Fuel is constantly fed from an external reservoir kept at 
a constant concentration. The fixed point in ( |35| ) is stable therefore small 
perturbations in the temperature field around the zero will tend to die out 
and return back to the zero, but what if the perturbation was big enough 
(a match is thrown into the fuel) to cause the term pcp\ t° De significant? 

Here the auto-catalytic nature of (|34|) would cause temperature and fuel 
consumption to increase inside the ignited area and the fire starts spreading 
across space, leaving regions with low concentrations of fuel. Given the feed 
parameters F and k the behavior is controlled by the value of Da- Large 
values of Da would cause fast moving fire wavefronts and small values of 
Da cause stable standing spots of fire fueled by the fast moving C. One of 



the more interesting patterns arising from (34) is that of replicating spots 
predicted in simulations j0|] and confirmed in experiment |l6|]. Given the 
right parameter values an initial large perturbation breaks up into standing 
spots, the fuel concentration at the middle of the spot decreases and causes 
the spot to divide into two smaller expanding replicas and so on. 



Analytic spot-like solutions to (^J) were found in 25] but only in the 
one-dimensional case and in the Da << 1 limit. 

The different two-dimensional patterns resulting from the Euler integra- 
tion of (^) were classified in [15]. 



5.2 Reproducing the replicating spots phenomenon in mi- 
croscopic simulations. 



The following equations : 

2 • 10" 5 V 2 pc - POP a + 0-018(1 - pc) 



d PC _ ln -:, r 2, 



dt 
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dp A _ , ^2 , , , .2 



at 



W- b V'pA + pc Pa ~ 0.074 PA (37) 



fall in the parameter space domain described by Pearson[15] to produce 
the replicating spots behavior. Pearson has used lattice-spacing of 0.01 and 
has started his simulation in the trivial state: pa = 0, pc = apart from 
the perturbated square put at: pa = 0.25, pc = 0.5. The reaction dynamics 
corresponding to equation (p7|) are: 



1. A + A + C^ A + A + A(l) 

2. A -» 0(0.74) 

3. 0(0.018) 

4. -> C(0.018) 

The simulation space (of size 64x64 sites) is seeded with 1000 Cs and 
no As, except for the perturbated square (of size 20x20 sites) that is seeded 
with 250 As per cell and 500 Cs per cell, no was chosen to be 10 7 and lattice 
spacing was chosen to be 0.01 (so a concentration of 1 corresponds to 1000 
reactants per lattice site). 

Although we have been successful in reproducing the general replicating 
spots behavior, a closer look at the results as shown in figures @ and (||) 
shows that the randomness introduced in the microscopic simulation has the 
effect of breaking the square symmetry preserved by the PDE simulations. 
Pearson originally achieved this symmetry breaking by introducing noise 
in the initial conditions, but the same effect is created by the microscopic 
simulation. Figure ( fTo) ) shows the results of the simulation at a later time, 
replication has formed more spots. Figure (|i~l| ) shows the results of a PDE 
simulation of the same system. We can see the similarity in the results of 
the two simulations, due to the fact that uq is large. We shall see, in the 
next subsection, that when other parameters are chosen (including no), the 
two systems give crucially different results. 

5.3 The uniqueness of persistently dynamic reaction-fronts 
to microscopic simulation. 

We now present an example in which microscopic simulations of the Gray- 
Scott model create different results than the ones obtained by the PDE 
approach. The phenomenon we are interested in is the presence persistent 
reaction fronts propagating through the simulation space. Much research 
has been devoted to reaction fronts ||§[14][24][12|. Patterns which are 
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Figure 8: Replicating spots of A created by the microscopic simulation. 
Spots in the process of dividing are shown. 
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Figure 9: The same simulation as in fig(||) at a later time, the original 
spots have split in two. 




Figure 10: The same simulation as in fig© and @ but at a later time. 
This figure shows species C. 
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Figure 11: This figure shows results of PDE simulation of the replication 
spots. Again showing species C. 

non-stationary in time and inhomogeneous in space at the same time are 
rare in a homogeneous medium. We suggest a system where the reaction- 
fronts are spatially and temporally non-stationary within a homogeneous 
medium (all the participating reactants have non-zero equal diffusion rates, 
so the system is not equivalent to a non-homogeneous medium system). 
Our persistent reaction-fronts are created by the following mechanism: A 
localized small A area consumes C reactants in its surroundings and creates 
reaction-fronts of high A concentration propagating across the simulation 
space cleaning areas from the presence of Cs. The reaction- front then runs 
out of high C concentration areas and decays, giving rise to the renewal of 
Cs and new wavefronts produced by the surviving As and so on. We have 
simulated microscopically the following reaction-diffusion system: 

1. A + A + C^ A + A + A(8) 

2. A -> 0(0.3) 

3. C -> 0(0.02) 

4. -> C(0.1) 
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Figure 12: Clear high A concentration reaction-fronts can be seen during 
microscopic simulations. These reaction-fronts are persistent, being sup- 
ported by surviving A reactants. 



The diffusion coefficients are Da = 1 and Dq = 1. The simulation size 
chosen was 80x80. uq = 1 and Ax = 1. The simulation space was seeded 
with 1 C per site, and no A's were seeded except for a perturbed square of 
size 64x64 sites where 5 reactants of species A per site were seeded. The 
simulation shows persistent reaction-fronts of the type seen in fig(12). 

Using Euler Integration to simulate the system we get different results, 
the initial square perturbation does propagate in the form of a reaction- 
front shown in fig(|l3|), but after the initial front cleans the space from high 
C concentration areas the front decays and gives way to the total elimination 
of As from the simulation space. Unlike the microscopic case in which the 
discretization and randomness leave behind some islands of A from which 
the next reaction-front can emerge, the Euler Integration drives the system 
towards the pa = steady-state. Introducing noise in the initial conditions 
does not change the situation the system tends to smooth out these noises 
and is driven to the constant steady state. 

We have tried to reproduce the persistent reaction-fronts pattern using 
Euler Integration with no success, although it is plausible that the same 
behavior will appear for high values of no- 



27 




Figure 13: The initial reaction-front seen when integrating the correspond- 
ing PDE system with Ax = 1 and no = 1. The reaction- front dies out and 
the system is driven to the constant pa = steady-state. 
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6 Microscopic simulation of Marketing models 



6.1 Using rate equations to describe the marketing of prod- 
ucts 

The application of physical sciences methods to economic and financial re- 



search, nowadays common practice among the scientific community! 17], was 



initiated by the work of Bachelier[ll] and Mandelbrot [|J] . The basic situa- 



tion is similar to the other areas of research we have previously discussed, 
global "macroscopic" economic phenomena are generated by the underly- 
ing "microscopic" process of buy and sell. Under these circumstances it is 
natural to try and use microscopic simulation of market models in order to 
reproduce as an example we could take the spreading of steam engines or 
gunpowder starting from localized innovative centers and sweeping across 
wide regions of the globe. 

When using microscopic simulation to describe product marketing one 
can regard a lattice-space element in two different ways: 

1. As a geographical region - neighborhood, town or country. 

2. As a business entity - company or corporate. 

In the first case neighboring elements represent geographically close regions, 
in the latter case they represent companies that are in business contact. 
The reactants could stand for any kind of valuable passing from one place 
or business entity to another. One could have reactants representing money, 
products, ideas, manpower or technology. Reactions could represent eco- 
nomic processes in which valuables are transformed, lost or created. Diffu- 
sion stands for the transfer of these valuables from one location or business 
entity to another. Although the use of microscopic simulations in the con- 
text of marketing seems natural enough there still are some flaws to point 
at and points to defend. The microscopic simulation is probabilistic in na- 
ture while some of the processes we try to represent are deterministic. As 
we explained in previous sections rates represent the occurrence rate of an 
underlying Poisson process for some of the economic processes described we 
have no reason to think that they are Poissonic in nature. Furthermore the 
microscopic simulation's lattice-space elements are always connected with 
exactly 4 neighboring sites, thus disabling the simulation of environments in 
which the degree of "connectedness" varies from site to site. If we consider 
sites to be business entities as we have previously suggested, we certainly 
would like to consider different amounts of connectivity between business 
entities - a feature not supported by the microscopic simulation system. 
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6.2 The "Tamagotchi" model 

In collaboration with J. Goldenberg and D. MazurskyQwe have devised the 
following set of reactions :^ 

1. B + C -> A + B + B(8) 

2. A — > 0(0.7) 

3. B — » 0(1.5) 

4. C(0.1) 

5. C -> 0(0.02) 

The inspiration to these reaction dynamics comes from the wave of "Tam- 
agotchi" games we have been subjected to during a short period of 1996. 
These reaction dynamics were simulated in a system of size 128x128 sites 
seeded with 1 reactant of species C only. A square of size 64x64 sites was 
seeded with 1 reactant of species C per cell and 5 reactant of species B per 
cite. Other coefficients are: no = 1, Ax = 1. 

• A - represents a product ("Tamagotchi"). 

• B - represents the idea or concept of the product. 

• C - represents money. 

In this model neighboring locations should be interpreted in the geograph- 
ical sense. The only diffusing reactant is B with diffusion coefficient 1 but 
since B represents an idea or concept the diffusion it performs is replica- 
tive. By replicative we mean that instead of hopping from one cell to a 
neighboring cell a copy of the original reactant is created and planted in 
the neighboring cell. This replicative diffusion represents the fact that ideas 
or concepts do not physically pass from one location to another, instead a 
location exposed to a new idea or concept in a neighboring location creates 
its own copy of the original idea. In terms of the microscopic simulation al- 
gorithm the only change is that when executing the diffusion of B, a B unit 
is added to the target location but none is subtracted in the original loca- 
tion. The first reaction represents the following process: A person exposed 
to the product concept that is in possession of money spends the money 

5 Hebrew University School of Business Management 

6 Numbers at the right of reactions are the corresponding reaction-rates. 
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Figure 14: Snap shot of A's concentration clear wavefronts can be seen. 
Sites with concentration above 2 are drawn black. 

on buying the product and a new concept or idea of the product is "born" 
in that persons mind. The second reaction represents the loss of products 
due to old age. The third reaction represents the fact that ideas tend to be 
forgotten. The fourth reaction represents an influx of money and the fifth 
reaction represents the spending of money on products other than A. The 
"replicative" diffusion that B undergoes represents the influence ideas have 
on neighboring locations. 

6.3 Wave fronts the "Tamagotchi" model 

During the simulations of the "Tamagotchi" model we have observed long 
lasting wavefronts of high A concentrations moving across the simulation 
space. The process giving rise to such waves seems to be the following: 
Starting with a small amount of A in a C rich area the concepts created 
by this small A concentration propagate to neighboring sites and induce 
the decrease of C concentration and increase in A concentration in these 
sites. The areas to which this B influence arrives then become C poor areas 
stopping the increase in A's concentration giving rise to their decrease until 
C's concentration is high enough to sustain the next wavefront. 
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Figure 15: C(r,to,t) plotted as a function of r for £=3,4,5, the higher the 
plot the lower is t. A clear 'hump' can be seen to advance in r as t increases, 
corresponding to the evolution of the wavefront. 



In order to estimate the propagation speed of the wavefronts we have 
used the following correlation integral. Let f2 denote our two-dimensional 
reaction chamber and let p^(x,t) denote the concentration of A at location 
x and time to for x € O and to > 0. For r > let: 

V(n) = [ ldx (38) 
Jn 



1 f 2n f ( rcosd \ 

V{9) 2 J n P A ^,to)dx J^p A (y,t + t)dy (39) 

Given t > and r > (|39|) is a measure for the correlation between A's 
concentration at time to and A's concentration at time to + 1 at locations 
with distance r. If pa is static C(r,to,t) is maximal for r = for all ts. If 
on the other hand pa is moving at speed v we expect for a given t that C 
is maximal for r=vt. 
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Figure 16: C(r, to + 50, t) plotted as a f unction of r for £=3,4,5, the higher 
the plot the lower is t. The same 'hump' as in (|l~5|) although to has increased 
by 50. 



In both (15) and fllq ) we can see an increase of 2 in r as t increases by 
2 giving rise to an estimate of: 



1 



(40) 



When speeds are measured in units of lattice-sites per simulation-time. 
The similar approximations due to the different choices of to indicate the 
constant speed of advance of the wavefront. 
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